In this code we'll optimize code used in #011, try to peform a greater number of simulations using log returns in a continuous function, and then convert the results only to arithimetical returns. Expected to reduce considerably time spending in monte carlo runs.
simplified calcs using risk free = 0, makes sharpe = rm/std
#import Libraries
import pandas as pd
from pandas_datareader import data as pdr
import numpy as np
from scipy.optimize import minimize
import random
import plotly.express as px
import plotly.graph_objects as go
import matplotlib.pyplot as plt
import matplotlib.ticker as mtick
def rand_weights(n):
''' Produces n random weights that sum to 1 '''
k = np.random.rand(n)
return k / sum(k)
def run_portfolio(df, weights, initial_investment, risk_free):
''' Performs an asset allocation and calc Portfolio Metrics'''
# The function returns:
# (0) Portfolio Data Frame
# (1) Expected portfolio return
# (2) Expected volatility
# (3) Sharpe ratio
# (4) Return on investment
# (5) Final portfolio value in dollars
# Perform asset allocation using the random weights (sent as arguments to the function)
portfolio_df = df.copy()
# Scale stock prices using the "price_scaling"
scaled_df = df.copy()
for i in df.columns[0:]:
scaled_df[i] = df[i]/df[i][0]
#enumerate method links Stocks tickers in columns along with a counter position weight (i), like an index
for i, stock in enumerate(scaled_df):
portfolio_df[stock] = weights[i] * scaled_df[stock] * initial_investment
# Sum up all values and place the result in a new column titled "portfolio value [$]"
portfolio_df['Total Value [$]'] = portfolio_df.sum(axis = 1, numeric_only = True)
# Calculate the portfolio percentage daily return and replace NaNs with zeros
portfolio_df['Daily Return [%]'] = portfolio_df['Total Value [$]'].pct_change(1) * 100
portfolio_df.replace(np.nan, 0, inplace = True)
# Calculate the return on the investment = last final value of the portfolio compared to its initial value
return_on_investment = ((portfolio_df['Total Value [$]'][-1:] - portfolio_df['Total Value [$]'][0])
/ portfolio_df['Total Value [$]'][0]) * 100
# Daily change of every stock in the portfolio
portfolio_daily_return_df = portfolio_df.drop(columns = ['Total Value [$]', 'Daily Return [%]'])
portfolio_daily_return_df = portfolio_daily_return_df.pct_change(1)
portfolio_daily_return_df.replace(np.nan, 0, inplace = True)
# Portfolio Expected Return formula
expected_portfolio_return = np.sum(weights * portfolio_daily_return_df.mean() ) * 252
# Portfolio standard deviation, volatility (risk)
covariance = portfolio_daily_return_df.cov() * 252
expected_volatility = np.sqrt(np.dot(weights.T, np.dot(covariance, weights)))
# Calculate Sharpe ratio
sharpe_ratio = (expected_portfolio_return - risk_free)/expected_volatility
return expected_portfolio_return, expected_volatility, sharpe_ratio, portfolio_df['Total Value [$]'][-1:].values[0], return_on_investment.values[0]
# file_name = input('Input the CSV file name: ')
# initial_investment = int(input('Input the initial investment: '))
# risk_free = float(input('Input the risk free annual rate: ')) # https://ycharts.com/ind icators/10_year_treasury_rate
# sim_runs = int(input('Input how many runs to perform on Monte Carlo Simulation: '))
file_name = 'BRL5' # CSV file with 5 Brazilian stocks and 10 year Ajusted Close Prices
initial_investment = 1000000 # 1mm
risk_free = 0.039 #3.9%
runs = 100000 # 100k
# read CSV file
df = pd.read_csv(file_name)
# replace the colum Date to Index
df.set_index(['Date'], inplace = True)
tickers = list(df.columns)
#Log is continuous and approx arithmetic daily return (252. over basis)
log_returns = df.pct_change().apply(lambda x: np.log(1+x)).dropna() #return for each stocks position column
mean_log_returns = log_returns.mean() # mean of Log returns | Mu
covariance_matrix = log_returns.cov() # covariance matrix | Sigma
#create arrays to store results
log_expected_return_array = np.zeros(runs) # expected return | E(R)
volatility_array = np.zeros(runs) # variances matrix | Sigma
sharpe_array = np.zeros(runs) #
weights_array = np.zeros((runs, len(df.columns))) # | W
# loop to perform monte carlo with random weights
for i in range(runs):
w = np.random.rand(len(df.columns))
w = w/np.sum(w)
weights_array[i, :] = w
log_expected_return_array[i] = np.sum(mean_log_returns * w * 252)
volatility_array[i] = np.sqrt(np.dot(w.T, np.dot(covariance_matrix*252, w)))
sharpe_array[i] = log_expected_return_array[i]/volatility_array[i]
# maximum sharpe ratio over simulation runs
opt_sharpe_index = sharpe_array.argmax()
opt_sharpe = sharpe_array.max()
opt_weights = weights_array[opt_sharpe_index]
# convert log return to arithmetic return
expected_return_array = np.exp(log_expected_return_array) - 1
#run for optimal portfolio
opt_return = np.exp(np.sum(mean_log_returns * opt_weights * 252)) - 1
opt_volatility = np.sqrt(np.dot(opt_weights.T, np.dot(covariance_matrix*252, opt_weights)))
# for y axis, select the maximum and minumum return points to plot efficient frontier
frontier_Y = np.linspace(expected_return_array.min(), expected_return_array.max(), 50)
# select each allocation, calculate and covert log to arithmetic return
def select_returns(w_1):
return_w1 = np.exp(np.sum(mean_log_returns * np.array(w_1)) * 252) - 1
return return_w1
# check if all capital is allocated into the assets, sum shoud be equal to 1
def check_allocation(w_1):
return np.sum(w_1)-1
# select each allocation, calculate volatility with w_1 and covariance matrix
def select_volatility(w_1):
w_1 = np.array(w_1)
w_1_volatility = np.sqrt(np.dot(w_1.T, np.dot(covariance_matrix*252, w_1)))
return w_1_volatility
# supose a equal weighted vector to start plot
w_0 = [1/len(tickers)] * len(tickers)
# impose bounds to plot function
bounds = tuple([(0, 1) for stock in tickers])
# set a vector to plot X axis
frontier_X = []
#run for max e min returns to plot X axis. verify if w sums to 1, if return selected appox to monte carlo return R
for R in frontier_Y:
# impose constraints to run least squares
constraints = ({'type':'eq', 'fun':check_allocation}, {'type':'eq', 'fun': lambda w: select_returns(w) - R})
# Minimize a scalar function of one or more variables using Sequential Least Squares Programming (SLSQP).
result = minimize(select_volatility,w_0,method='SLSQP', bounds = bounds, constraints = constraints)
frontier_X.append(result['fun'])
# Create a DataFrame that contains volatility, return, and Sharpe ratio for all simualation runs
sim_out_df = pd.DataFrame({'Volatility': volatility_array.tolist(), 'Portfolio_Return': expected_return_array.tolist(), 'Sharpe_Ratio': sharpe_array.tolist() })
# Plotting data
fig1 = px.line(df, x= frontier_X, y= frontier_Y)
fig = px.scatter(sim_out_df, x = 'Volatility', y = 'Portfolio_Return', color = 'Sharpe_Ratio', size = 'Sharpe_Ratio', hover_data = ['Sharpe_Ratio'], size_max=10 )
fig.add_trace(go.Scatter(x = [opt_volatility], y = [opt_return], mode = 'markers', name = 'Optimal Point', marker = dict(size=[20], color = 'green')))
fig2 = go.Figure(data=fig.data + fig1.data)
fig2.update_layout(coloraxis_colorbar = dict(y = 0.7, dtick = 3),
# title = 'porfolio results ',
xaxis_title = 'Volatility(STD)',
yaxis_title = 'Return',
xaxis=dict(tickformat=".1%"),
yaxis=dict(tickformat=".1%"))
fig2.update_layout({'plot_bgcolor': "white"})
fig2.show()
print('Best Portfolio Metrics Based on {} Monte Carlo Simulation Runs:'.format(runs))
print(" - Portfolio Stocks = {}".format(tickers))
print(" - Portfolio Weights = {} %".format(opt_weights.round(4)*100))
print(' - Portfolio Expected Annual Return = {:.02f} %'.format(opt_return * 100))
print(' - Portfolio Standard Deviation (Volatility) = {:.02f} %'.format(opt_volatility * 100))
print(' - Sharpe Ratio = {:.02f}'.format(opt_sharpe))